MATLAB之LU分解法（十）

您所在的位置：网站首页 › matlab 对焦 › MATLAB之LU分解法（十）

MATLAB之LU分解法（十）

2024-07-06 12:47| 来源: 网络整理| 查看: 265

LU分解 1.LU分解的基础知识

矩阵的LU分解又称为矩阵的三角分解，即将一个矩阵分解为一个下三角矩阵L和一个上三角矩阵U，即 A = L U A=LU A=LU,其在方程组的求解和求矩阵的逆有许多应用。

LU分解的求解命令是lu，基本使用格式如下：

调用格式说明[L,U]=lu(A)将一个矩阵分解为一个下三角L和一个上三角U[L,U,P]=lu(A)对A进行分解，其中L为单位下三角矩阵，U为上三角矩阵，P为置换矩阵，满足LU=PA

例1,对矩阵 A = [ 1 2 3 4 5 6 7 8 2 3 4 1 7 8 5 6 ] A=\left[\begin{array}{llll}1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ 2 & 3 & 4 & 1 \\ 7 & 8 & 5 & 6\end{array}\right] A=⎣⎢⎢⎡1527263837454816⎦⎥⎥⎤进行分解。

A=[1,2,3,4;5,6,7,8;2,3,4,1;7,8,5,6]; [L1,U1]=lu(A) L1 = 0.1429 1.0000 0 0 0.7143 0.3333 1.0000 0 0.2857 0.8333 0.2500 1.0000 1.0000 0 0 0 U1 = 7.0000 8.0000 5.0000 6.0000 0 0.8571 2.2857 3.1429 0 0 2.6667 2.6667 0 0 0 -4.0000 [L2,U2,P]=lu(A) L2 = 1.0000 0 0 0 0.1429 1.0000 0 0 0.7143 0.3333 1.0000 0 0.2857 0.8333 0.2500 1.0000 U2 = 7.0000 8.0000 5.0000 6.0000 0 0.8571 2.2857 3.1429 0 0 2.6667 2.6667 0 0 0 -4.0000 P = 0 0 0 1 1 0 0 0 0 1 0 0 0 0 1 0

从上面两者的计算结果比较分析可以得出：使用第一种方法得到的上三角矩阵并不是严格的上三角矩阵，这样对于以后的分析和计算显然是不利的，所以在工程计算或者其他计算中都采用第二种方法。

2.LU分解法求解方程组

LU分解法求解的基本思路就是：将系数矩阵A进行分解，得到 L U = P A LU=PA LU=PA,然后利用 L y = P b Ly=Pb Ly=Pb,最后再解 U x = y Ux=y Ux=y,从而求解方程组的解。

例2：利用LU分解法求解方程组 { x 1 + x 2 − 3 x 3 − x 4 = 1 3 x 1 − x 2 − 3 x 3 + 4 x 4 = 4 x 1 + 5 x 2 − 9 x 3 − 8 x 4 = 0 \left\{\begin{array}{l}x_{1}+x_{2}-3 x_{3}-x_{4}=1 \\ 3 x_{1}-x_{2}-3 x_{3}+4 x_{4}=4 \\ x_{1}+5 x_{2}-9 x_{3}-8 x_{4}=0\end{array}\right. ⎩⎨⎧x1+x2−3x3−x4=13x1−x2−3x3+4x4=4x1+5x2−9x3−8x4=0 的解。

A = [1,1,-3,-1;3,-1,-3,4;1,5,-9,-8]; b=[1;4;0] b = 1 4 0 [L,U,P]=lu(A) L = 1.0000 0 0 0.3333 1.0000 0 0.3333 0.2500 1.0000 U = 3.0000 -1.0000 -3.0000 4.0000 0 5.3333 -8.0000 -9.3333 0 0 0 0.0000 P = 0 1 0 0 0 1 1 0 0 y = L'*P*b y = 4.3333 0.2500 1.0000 >> x = U'*y x = 13.0000 -3.0000 -15.0000 15.0000

利用函数实现如下：

function x= solvebyLU(A,b) %此函数用于利用LU分解法解方程组 flag=isexist(A,b); %判断方程组是否有解 if flag==0 disp('方程组无解'); x = []; return; else r = rank(A); [m,n]=size(A); [L,U,P]=lu(A); b = P*b; %解Ly=b y(1) = b(1); if m>1 for i=2:m y(i)=b(i)-L(i,1:i-1)*y(1:i-1)'; end end y = y'; %解Ux=y得方程组的解 x0(r)=y(r)/U(r,r); if r>1 for i=r-1:-1:1 x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i); end end x0 = x0'; if flag==1 x=x0; return; else format rat; Z = null(A,'r'); [mz,nz]=size(Z); x0(r+1:n)=0; for i=1:nz t=sym(char([107 48+i])); k(i)=t; end x=x0; for i=1:nz x = x+k(i)*Z(:,i); end end end

运行结果：

A = [1,1,-3,-1;3,-1,-3,4;1,5,-9,-8]; b = [1;4;0]; x = solvebyLU(A,b) x = (3*k1)/2 - (3*k2)/4 + 5/4 (3*k1)/2 + (7*k2)/4 - 1/4 k1 k2

【本文地址】

MATLAB之LU分解法（十）

MATLAB之LU分解法（十）

今日新闻

推荐新闻